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Abstract 

We present an iterative algorithm for calculating approximate greatest common 
divisor (GCD) of univariate polynomials with the real or the complex coeffi- 
cients. For a given pair of polynomials and a degree, our algorithm finds a pair 
of polynomials which has a GCD of the given degree and whose coefficients 
are perturbed from those in the original inputs, making the perturbations as 
small as possible, along with the GCD. The problem of approximate GCD is 
transfered to a constrained minimization problem, then solved with the so-called 
modified Newton method, which is a generalization of the gradient-projection 
method, by searching the solution iteratively. We demonstrate that, in some 
test cases, our algorithm calculates approximate GCD with perturbations as 
small as those calculated by a method based on the structured total least norm 
(STLN) method and the UVGCD method, while our method runs significantly 
faster than theirs by approximately up to 30 or 10 times, respectively, compared 
with their implementation. We also show that our algorithm properly handles 
some ill-conditioned polynomials which have a GCD with small or large leading 
coefficient. 

Keywords: Approximate polynomial GCD, Gradient-projection method, 
Ill-conditioned problem, Optimization. 



1. Introduction 

For algebraic computations on polynomials and matrices, approximate alge- 
braic algorithms are attracting more attention than before. These algorithms 
take inputs with some "noise" such as polynomials with floating-point number 
coefficients with rounding errors, or more practical errors such as measurement 
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errors, then, with minimal changes on the inputs, seek a meaningful answer 
that reflect desired property of the input, such as a common factor of a given 
degree. By this characteristic, approximate algebraic algorithms are expected 
to be applicable to more wide range of problems, especially those to which exact 
algebraic algorithms were not applicable. 

As an approximate algebraic algorithm, we consider calculating the approx- 
imate greatest common divisor (GCD) of univariate polynomials with the real 
or the complex coefficients, such that, for a given pair of polynomials and a 
degree d, finding a pair of polynomials which has a GCD of degree d and whose 
coefficients are perturbations from those in the original inputs, while making the 
perturbations as small as possible, along with the GCD. This problem has been 
extensively studied with various approaches inclu ding the Euclidean method on 
the polynomial remainder sequence (PRS) f fill. |22[ 23 1 ) , the singular value 
decomposition (SVD) of the Sylvester matrix (p}, [10|). the LU or QR factor- 
ization of the Sylvester and/or Bezout matrix or their displacements 
l32ll V\ , Pade approximation ([H}), optimization strategies (Q, 
|31j). Furthermore, stable methods for ill-conditioned problems have 
been discussed (0, @, [13, [H|). 





Among methods in the above, we focus our attention on optimization strat- 
egy in this paper, especially iterative method for approaching an optimal so- 
lution, after transferring the approximate GCD problem into a constrained 
minimization problem. Already proposed algorithms utilize iterative methods 
including the Levenberg-Marquardt method (|6|), the Gauss-Newton method 
(Uni) and the structured total least norm (STLN) method ([HI], Q). Among 
them, STLN-based methods have shown good performance calculating approx- 
imate GCD with sufficiently small perturbations efficiently. 

Here, we utilize the so-called modified Newton method ([13]), which is a gen- 
eralization of the gradient-projection method ( 20] ) , for solving the constrained 
minimization problem. This method has interesting features such that it com- 
bines the projection and the restoration steps in the original gradient-projection 
method, which reduces the number of solving a linear system. We demonstrate 
that our algorithm calculates approximate GCD with perturbations as small as 
those calculated by the STLN-based methods, while our method show signifi- 
cantly better performance over them in its speed compared with their imple- 
mentation, by approximately up to 30 times. Furthermore, we also show that 
our algorithm can properly handle some ill-conditioned problems such as those 
with GCD containing small or large leading coefficient. We call our algorithm 
GPGCD after the initials of the gradient projection method. 

In this paper, we present the following expansion from the previous results 
((26[, [25|]) as presenting the new algorithm for monic polynomials to calculate 
perturbed polynomials without giving perturbations for the leading coefficients; 



1 Note that the article by Bini and Boito [2| has absence of reference to a literature on 
computation on structured matrices by Pan [lSl , whereas the dissertation by Boito Q has no 
such omission. 
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providing experiment results for new test polynomials that have been prepared 
more carefully and comparison with the UVGCD method ([3(|) (in Section [572)) ; 
adding more experiments for comparison of our algorithm with the STLN-based 
method and the UVGCD method (in Section [575)1 . 

The rest part of this chapter is organized as follows. In Section [5J we trans- 
form the approximate GCD problem into a constrained minimization problem. 
In Section [3) we review the framework of the gradient-projection method and 
the modified Newton method. In Section H) we show an algorithm for cal- 
culating the approximate GCD, and discuss issues in the application of the 
gradient-projection method or the modified Newton method. In Section [5) we 
demonstrate performance of our algorithm with experiments. 



2. Formulation of the Approximate GCD Problem 



Let F(x) and G(x) be univariate polynomials with the real or the complex 
coefficients, given as 



F(X) = f m X m + /, 

G(x) = g n x n + g n -\x 



\x 

ra-l 



+ /o, 

go, 



(i) 



with < n < m. We permit F and G to be relatively prime in general. For a 
given integer d satisfying < d < n, let us calculate a deformation of F(x) and 
G(x) in the form of 



F{x) = F(x) + AF(x) = H(x) ■ F(x), 
G(x) = G{x) + AG(x) = H(x) ■ G(x), 



(2) 



where AF(x), AG(x) are polynomials whose degrees do not exceed those of 
F(x) and G(x), respectively, H(x) is a polynomial of degree d, and F{x) and 
G(x) are pairwise relatively prime. If we find F, G, F, G and H satisfying 
([2]), then we call H an approximate GCD of F and G. For a given degree d, 
we tackle the problem of finding an approximate GCD H while minimizing the 
norm of the deformations || AF(x)||| + II^G(x)|||. 

To make the paper self-contained, we define notations in the theory of sub- 
resultants used below. 

Defnition 1 (Sylvester Matrix). Let F and G be defined as in fTJ). The Sylvester 
matrix of F and G, denoted by N(F,G), is an (m + n) x (to + n) matrix 
constructed from the coefficients of F and G, such that 



\ 



N(F,G) = 



go 



fo 



go J 
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Defnition 2 (Subresultant Matrix). Let F and G be defined as in ([lj. For 

< j < n, the j-th subresultant matrix of F and G, denoted by Nj(F, G), is an 
(m + n — j) x (m + n — 2j) sub-matrix of N(F, G) obtained by taking the left 
n — j columns of coefficients of F and the left m — j columns of coefficients of 
G, such that 



Nj(F,G) 



ffn 



fo 



fo 



9n 



9o 



\ 



9oJ 



(3) 



Defnition 3 (Subresultant). Let F and G be defined as in (fTJ). For < j < n 

and k = 0, . . . , j, let Nj^ = Nj t k(F, G) be a sub-matrix of Nj(F, G) obtained 
by taking the top m + n — 2j — 1 rows and the (m + n — j — fc)-th row (note 
that Nj y k{F,G) is a square matrix). Then, the polynomial 



S i (F,G) = |iV iJ |^- 
is called the j-th subresultant of F and G. 



Now, in the case F(x) and G(x) have a GCD of degree d, then the theory 
of subresultants tells us that the (d — l)-th subresultant of F and G becomes 
zero, namely we have 

S d _i(F,G) = 0. 

Then, the (d — l)-th subresultant matrix Nd-i(F, G) has a kernel of dimension 
equal to 1. Thus, there exist polynomials A{x), B(x) G R[x] or C[x] satisfying 



AF + BG = 0, 



(4) 



with deg(A) < n — d and deg(_B) < m — d and A(x) and B(x) are relatively 
prime. Therefore, for the given F(x), G{x) and d, our problem is to find AF(x), 
AG{x), A(x) and B{x) satisfying Eq. (gj while making ||Z\F|||+||Z\G||| as small 
as possible. 



2.1. The Real Coefficient Case 

Assuming that we have F(x) and G[x) as polynomials with the real co- 
efficients and find an approximate GCD with the real coefficients as well, we 
represent F(x), G(x), A(x) and B(x) with the real coefficients as 

F(x) = f m x m + ■■■ + / x°, G(x) = g n x n + ■■■+ g x°, 

A{x) = a n _ d x n - d + ■■■+ a x°, B(x) = b m - d x m - d + ■■■ + b x°, 
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respectively, thus HAFHf. + ||ZiC||| and Eq. (0} become as 

||AF||| + \\AGf 2 = (f m - f m f + ... + (/„_ / ) 2 + (g n - g n f + ... + (g - g )\ 

(6) 

N d ^(F,G)-v = 0, (7) 

respectively, with Nj(F,G) as in and 

v = *(a„_ d , ...,a , b m -d, ■ ■ ■ ,b )- (8) 

Then, Eq. ([7]) is regarded as a system ofm + n — d + 1 equations in f m , . . . , /o, 
9m ■ ■ ■ , go, a n-d, ■ ■ ■ j ao, b m -d, . . . ,bo, as 

<Zl = fm a n-d + 9nb m -d = 0, • • ■ , Om+n-d+1 = /o a o + ffO&O = 0, (9) 

by putting as the j-th row. Furthermore, for solving the problem below 
stably, we add another constraint enforcing the coefficients of A(x) and B(x) 
such that + \\B(x)\\l = 1; thus we add 

9o = 4-d + • • • + ao + &m-d + • • • + b 2 - 1 = (10) 

into Eq. ©. 

Now, we substitute the variables 

{fm, ■ ■ ■ , fo,9n, ■ ■ ■ ,9o, a n -d, • ■ • , ao, b m -d, ■ ■ ■ ,bo) (11) 
as x — (xi, . . . ,X2tm+n-d+2)), thus Eq. © and §§§ with (|TU)) become 

f(x) = (si - / m ) 2 H h (x m+1 - /o) 2 

+ (x m +2 - 9n) 2 H h (a; m+ „ +2 - go) 2 , (12) 

= *(go(^), 9i ■ • ■ , gm+n-d+it*)) = 0, (13) 

respectively. Therefore, the problem of finding an approximate GCD can be 
formulated as a constrained minimization problem of finding a minimizer of the 
objective function f(x) in (fT2|) . subject to q(x) = in Eq. fp~3|) . 

2.2. The Complex Coefficient Case 

Now let us assume that we have F(x) and G(x) with the complex coefficients 
in general, represented as 

F(x) = (f mA + f m ,2i)x m + ■■■ + (fo.l + fo, 2 i), 
G(x) = (g n ,i + g n ai)x n + H h (#o,i + go,2i), 

where /,• i, 2, 9j, 2 are real numbers; 1, and gj 1 represent the real parts; 
fj 2, 5j,2 represent the imaginary parts, with i as the imaginary unit, and find 
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an approximate GCD with the complex coefficients. Then, we represent F(x), 
G(x), A(x) and B(x) with the complex coefficients as 



(14) 



F(x) = (f mA + f m , 2 i)x m + ■■■ + (,/o,i + f , 2 i)x°, 

G(x) = (g n> i + g n ,2i)x n H h {g x° + g , 2 i)x°, 

A{x) = (a n - d ,i + a n - d , 2 i)x n ~ d H h (a ,i + a , 2 i)x° , 

B(x) = (& m _ d ,i + b m ^ da i)x m ~ d + ■■■ + (6 ,i + b 0j2 i)x°, 

respectively, where fj t %, fj t2 , §3,2, aj,2, bj t2 are real numbers. 
For the objective function, ||AF||2 + || ^C|l | becomes as 



£[(/,,! - fj.lf + (fj,2 h2) 2 } + - 9j,l? + - .9„2) 2 ]. (15) 

i=o j=o 
For the constraint, Eq. Q becomes as 

/7m,l + /m,2« <7n,l + <?n,2« 



/o,l + /o,2« /m,l + /m,2» <70,1 + #0,2* <?n,l + <7ra,2* 



/o.l + /o,2* 



.90,1 +5o,2V 



«o,i + ao,2« 

u rn — a.l 

+ b 



\ &o,i+&o,2« J 



= 0. (16) 



By expressing the subresultant matrix and the column vector in (fl"6|) separated 
into the real and the complex parts, respectively, we express (1161) as 



with 



JVi 



m,l 



{N x +N 2 i)(v 1 +v 2 i) = 0, 

Jn,l \ ( fm,2 



(17) 



fo,l fm,l 90,1 9n,l 



, N 2 



9n,2 



fo,2 fm,2 §0,2 g n ,2 



/o,l ffO,l/ V /o,2 

t>l = t {a n -d,l, ■ ■ ■ , <JQ,1> b m -d,l, ■ ■ • , &0,l)> 
1>2 = t (a n -d,2, ■ ■ ■ , 0,Q,2t b rn -d,2, ■ ■ ■ , bo, 2 )- 



(18) 



50,2/ 
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We can expand the left-hand-side of Eq. (fT?) as 

(Ni + N 2 i){v 1 + v 2 i) = (NxVi - N 2 v 2 ) + i(N x v 2 + N 2 Vi), 
thus, Eq. (TIT)) is equivalent to a system of equations 



N 2 v 2 


= 0, 


Niv 2 - 


f N 2 v x = 0, 


(N x 


~N 2 ^ 


)(::) 


= 0. 


\N 2 


Ni , 





which is expressed as 

/ AT. _ AT~\ /-.,. \ 

(19) 

Furthermore, as well as in the real coefficients case, we add another con- 
straint for the coefficient of A(x) and B(x) as 

\\A(x)\\l + \\B(x)\\l = (4_ dl + • • • + a 2 0tl ) + (b 2 m _ d<1 + ■■■ + b 2 0A ) 

+ (al-d,2 + ■■■ + < 2 ) + (b 2 m - d ,2 + ■■■ + b 2 , 2 ) - 1 = 0, (20) 
which can be expressed together with (fT§|) as 

' l v x l v 2 -1\ f Vl \ 

JVi -N 2 « 2 = 0, (21) 
N 2 N x / \ 1 J 

where Eq. ([2H)l has been put on the top of Eq. (ITO)) . Note that, in Eq. (f2"Tj) , we 

have total of 2(m + n — d + 1) + 1 equations in the coefficients of polynomials 
in (|14[) as a constraint, with the j-th row of which is expressed as qj = 0, as 
similarly as in the real case © with (TTTTj) . 

Now, as in the real case, we substitute the variables 

(/m,l, • ■ • , fo.l, 9n,l, ■ ■ ■ j 90,1, fm,2, • • ■ , /o,2i <?n,2, ■ • ■ ,§0,2, 

Qln-d,l; • • • ) a 0,l, b m -d,l, ■ ■ ■ , &0,1) 0>n-d,2i ■ ■ ■ , a 0,2, &m— d,2; ■ • ■ , &0,2) (22) 

as x — (xt, ■ ■ ■ , ^4(m+n-d+2))) thus Eq. (TTS|) and (j2"Tj) become as 

f(x) =( Xl - fm.x) 2 + ■■■ + (X m+1 - / ,i) 2 

+ {Xm+2 - gn.l) 2 H h (£m+n+2 ~ 50,l) 2 

+ (^m+n+3 — /m,2) 2 + " " " + (x 2 m+n+3 ~ fo,2) 2 

) 2 , (23) 

g(cc) = *(gi(cc), . . . , fe(m+n-d+l)+l( aj )) = °) ( 24 ) 

respectively. Therefore, the problem of finding an approximate GCD can be 
formulated as a constrained minimization problem of finding a minimizer of the 
objective function f(x) in Eq. (|23l) . subject to q(x) — in Eq. 
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3. The Gradient-Projection Method and the Modified Newton Method 



In this section, we consider the problem of minimizing an objective function 
f(x) : R™ — > R, subject to the constraints q{x) = for q(x) = t (q 1 (x), q%{x), . . . , q m (x)), 
with m < n, where Qj(x) is a function of R™ — > R, and f(x) and qj(x) are twice 
continuously differentiable (here, we refer presentations of the problem to Tan- 



abe [24| and the references therein). 



If we assume that the Jacobian matrix 

J q {x) 

is of full rank, or 



dxj 



rank(Jq(a;)) = to, (25) 
on the feasible region V q defined by 

V q = {x e R" | q(x) = 0}, 

then the feasible region V q is an (n — TO,)-dimensional differential manifold in 
R" and / is differentiable function on the manifold V q . Thus, our problem is 
to find a point in V q , which will be a candidate of a local minimizer, satisfying 
the well-known "first-order necessary conditions" (for the proof, refer to the 
literature on optimization such as Nocedal and Wright [H}). 

Theorem 1 (First-order necessary conditions). Suppose that x* 6 V q is a 
local solution of the problem in the above, that the functions f(x) and q(x) are 
continuously differentiable at x* , and that we have (|25p at x* . Then, there exist 
a Lagrange multiplier vector A* € R m satisfying 

Vf{x*)-\j q {x*))\* = 0, q(x*) = 0. □ 

3.1. The Gradient- Projection Method 

Let Xk & R" be a feasible point, or a point satisfying Xk £ V q . Rosen's 
gradient projection method ([20() is based on projecting the steepest descent 
direction onto the tangent space of the manifold V q at Xk, which is denoted to 
T Xk and represented by the kernel of the Jacobian matrix J q (xk) as 

T Xk = ker(J q (x k )) = {ze R" | J q (x k )z = e R m }. (26) 

We have steepest descent direction of the objective function / at x k as 

Then, the search direction d k is defined by the projection of the steepest descent 
direction of / in ([27]) onto T Xk in ([26]) as 



d k = -P(x k )Vf(x k ). (28) 
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Here, P(xk) is the orthogonal projection operator on T Xk defined as 

P(x k ) = I - (J q (x k )) + (J q (x k )), 

where / is the identity matrix and (J q (xk)) + is the Moore-Penrose inverse of 
(Jq(xk))- Under the assumption (J5S]), we have 

(J q {Xk)) + = \jq(x k )) ■ (Jq(Xk) ■ HjgiXk)))' 1 

(see Tanabe O, Eq. (8)]). 

With an appropriate step width a k (see Remark [1]) satisfying < a k < 1, 

let 

Uk = x k + ak ■ d k . 

Since V q is nonlinear in general, yk may not in V q : in such a case, we take a 
restoration move to bring yk back to V q , as follows. Let x 6 R™ be an arbitrary 
point. Then, at yk, the constraint q{x) can be linearly approximated as 

q(y k + x) ~ q(y k ) + J q (y k )x. 

Assuming y k + x G V q , we have q(yk + x) = thus the approximation of x can 
be calculated as 

x = -{J q {Vk)) + q{yk). (29) 
If yk is sufficiently close to V q , then we can restore yk back onto by applying 
(|29|) iteratively for several times. Note that the restoration move can also be 
used in the case the initial point of the minimization process is away from the 
feasible region V q . 

Summarizing the above, we obtain an algorithm for the gradient projection 
as follows. 

Algorithm 1 (The gradient-projection method ([20])). 

Step 1 [Restoration] If the given point Xq does not satisfy Xq E V q , first 
move Xo onto V q by the iteration of Eq. (|29p , then let Xq be the restored 
point on V q . Let k = 0. 

Step 2 [Projection] For x k , calculate d k = -P(x k )V f{x k ) by (gSJ). If ||d fc || 
is sufficiently small for an appropriate norm, go to Step 4. Otherwise, 
calculate the step width ak by an appropriate line search method (see 
Remark [T]) then let yk,o = x k + afcdfc- 

Step 3 [Restoration] If q(yk,a) ^ 0, move y k) o back onto V q iteratively by 
([29]). Let y k j+i = y k ,i - (J q (yk,i)) + q(y k ,i) for I = 0, 1,2, . . .. When y Kl 
satisfies q{y kt i) ~ 0, then let x k+ i = yk,i and go to Step 2. 

Step 4 [Checking the first-order necessary conditions] If Xk satisfies 
Theorem [TJ then return Xk- 

Remark 1. Choosing appropriate step width in the iteration is a fundamental 
issue in optimization method and is discussed in standard literature of opti- 



mization {e.g. 16]). Although we simply set ak — 1 in our implementation, 



more sophisticated calculation of step width might improve accuracy and/or 
convergence of the algorithm (see also concluding remarks (Section [B])). 
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3.2. The Modified Newton Method 

The modified Newton method by Tanabe 24] is a generalization of the 
Newton's method, which derives several different methods, by modifying the 
Hessian of the Lagrange function. A generalization of the gradient-projection 
method combines the restoration step and the projection step in Algorithm [TJ 
For x k £ V q , we calculate the search direction d k , along with the associated 
Lagrange multipliers A fc+1 , by solving a linear system 



t 



(J q (x k ))\ ( d k \ _ [Vf(x k )\ 
O )\X k+ x)- \q(x k ))> {30 > 



then put x k +i = x k + a k ■ d k with an appropriate step width a k . Solving Eq. 
(|30|) under assumption (|25|) . we have 



d k = -P(x k )Vf(x k ) - (J q (x k )) + q(x k ), 
A fe+ i = t ((J (J (a ;fe )) + )V/( a;fc ) - (J q (x k ) ■ '(^(x*))) -1 ^**)- 

Note that, in d k in (pT|) . the term —P{x k )Vf(x k ) comes from the projection 
(|28|) . while another term — ( J q (x k )) + q(x k ) comes from the restoration ([29]) . 
If we have x k £ V q , the iteration formula (|30|) is equivalent to the projection 
(|28j) . After an iteration, the new estimate x k +i may not satisfy x k +\ £ V q : in 
such a case, in the next iteration, the point will be pulled back onto V q by the 
— (J q (xk)) + q(x k ) term. Therefore, by solving Eq. (|5U| iteratively, we expect 
that the approximations x k moves toward descending direction of / along with 
tracing the feasible set V q . 

Summarizing the above, we obtain an algorithm as follows. 

Algorithm 2 (The modified Newton method ([Hj])). 

Step 1 [Finding a search direction] For x k , calculate d k by solving the 
linear system ([30]) . If is sufficiently small, go to Step 2. Otherwise, 
calculate the step width a k by an appropriate line search method (see 
Remark]!]), let x k+ i = x k + a k d k , then go to Step 1. 

Step 2 [Checking the first-order necessary conditions] If x k satisfies 
Theorem Q] with sufficient accuracy, then return x k . 



4. The Algorithm for Approximate GCD 

In applying the gradient-projection method or the modified Newton method 
to the approximate GCD problem, we discuss issues in the construction of the 
algorithm in detail, such as 

• Representation of the Jacobian matrix J q {x) (Section 14. ip . 

• Stability of the algorithm by certifying that J q (x) has full rank (Section 

S3), 
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• Setting the initial values (Section 

• Regarding the minimization problem as the minimum distance problem 
(Section [O}, 

• Calculating the actual GCD and correcting the coefficients of F and G 
(Section EH)]). 

as follows. After presenting the algorithm, we give a modification for preserving 
monicity for the real coefficient case and running time analysis, and end this 
section with examples. 



4-1. Representation of the Jacobian Matrix 

For a polynomial P(x) € R[x] or C[a;] represented as 

P(x) = Pn x n + ---+ Po x°, 

let Cfe(P) be a complex (n + k, k + 1) matrix defined as 

(Pn \ 



C k {P) 



PO Pr, 



Po/ 



fc+1 



We show the Jacobian matrix in the real and the complex coefficient cases, 
both of which can easily be constructed in every iteration in Algorithms Q] and 

m 



4-1.1. The Real Coefficient Case 

For co-factors A{x) and B{x) as in ([5]), consider matrices C m (A) and C n (B). 
Then, by the definition of the constraint (fT3|) . we have the Jacobian matrix J q (x) 
(with the original notation of variables for x as in ([TTj) ) as 

/ 2 • *v \ 

Jq{x) ~ \C m (A) C n {B) Nd-i(F,G)J ' (32) 

with Nj(F,G) as in ([3]) and v as in ((8|), respectively. Note that the matrix 
J q {x) has m + n — d + 2 rows and 2(m + n — d + 2) columns. 

4-1.2. The Complex Coefficient Case 

For co-factors A(x) and -B(a;) as in (fT4|) . consider matrices C rra ( J 4) and C n (B) 
and express them as the sum of matrices consisting of the real and the imaginary 
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parts of whose elements, respectively, as 

/ 0>n—d,l \ ( a n-d,2 



C m {A) 



ao.i 



In-dA 



V "0,1 / 

C m (A) 1 + iC m (A) 2 , 

(bm-d,l \ 



ao,2 



C n (B) 



bo, 



b m -d,i 



fbm-d' 



bo,: 



V b 0A j 

= C„( J B) 1 +iC n (B) 2 , 

respectively, and define 

/ 0>n-d,l 



A 1 = [C m (A)i C n (fl)i] = 



Jm-d,l 



\ 

( a-n-d,2 



00,1 



A 2 = [C m (A) 2 C n (B) 2 ] 



-d.2 



ao,2 dn-d,2 bo,: 



V 



ao,2 



0>n-d,2 



°0,2 / 



bm-d:. 



b ,2 j 



0-0,1 a>n-d,i bo,i b m -d,i 



bo,i j 
\ 



-d.2 



bo,2 j 



(33) 



(Note that A\ and A 2 are matrices of the real numbers of m + n — d + 1 rows 
and m + n + 2 columns.) Then, by the definition of the constraint (pM)) . we have 
the Jacobian matrix J q {x) (with the original notation of variables for x as in 
(|22])) as 

/ 2 • Vl 2 ■ t v 2 \ 
J q (x) = A 1 -A 2 Ni -N 2 , (34) 
\A 2 A l N 2 Ni J 

with Ai and A 2 as in (|33|) and iVi, A?2, «i and «2 as in (fT8|) , respectively. 



4-2. Stability of the Algorithm 

In this paper, we treat the notion of "stability" of the algorithm as to keep 
that the Jacobian J q {x) in Algorithms [1] and [2] has full rank, whereas we usu- 
ally discuss stability as a notion in backward and/or forward error analysis of 
numerical algorithms [Tlj |. 
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In executing Algorithm[TJor[2J we need the algorithm to be stable in the sense 
that we need to keep that J q (x) has full rank: otherwise, we cannot correctly 
calculate (J q (x)) + (in Algorithm [lj or the matrix in (|30l) becomes singular (in 
Algorithm [5]) thus we are unable to decide proper search direction. For this 
requirement, we have the following observations. 

Proposition 1. Let x* G V q be any feasible point satisfying Eq. (|13[) . Then, if 
the corresponding polynomials do not have a CCD whose degree exceeds d, then 
J q (x*) has full rank. 

Proof. We prove the proposition in the real and the complex coefficient cases 
separately. 

4-2.1. The Real Coefficient Case 

Let x* = (/,„, . . . , /o, g n , ■ ■ ■ ,9o, a>n-d ■ • • ,a , b m - d , ...,b ) with its polyno- 
mial representation expressed as in ([5]) (note that this assumption permits the 
polynomials F(x) and G(x) to be relatively prime in general). To verify our 
claim, we show that we have rank(J q (x*)) = m + n — d + 2 with J q (x*) as in 
(|3^|) . Let us express J g (x*) — (Jl | Jr), where Jl and Jr are column blocks 
expressed as 

JL = (cJU C n (B)) ' Jr = (iV rf _i(F, G)) ' 

respectively. Then, we have the following lemma. 
Lemma 1. We have rank(JL) = m + n — d + 1. 
Proof. Let us express Jl = ( Jll | Jlr) > where 

• /LL = (a^)) ' = {cn(B)) ' 

and let Jl be a submatrix of Jl by taking the right m — d columns of Jll and 
the right n — d columns of Jlr- Then, we see that the bottom m + n — 2d 
rows of Jl is equal to N(A, P>), the Sylvester matrix of A(x) and B{x). By the 
assumption, polynomials A(x) and B(x) are relatively prime, and there exist 
no nonzero elements in Jl except for the bottom m + n — 2d rows, we have 
rank(JL) = m + n — 2d. 

By the above structure of Jl and the lower triangular structure of Jll and 
Jlr, we can take the left d+1 columns of Jll or Jlr satisfying linear indepen- 
dence along with the m + n — 2d columns in Jl. Therefore, these m + n — d + 1 
columns generate a (m+n — d+l)-dimensional subspace in pj"+ n - d + 2 satisfying 

f(n,. ■ • , x m+n - d+2 ) g K m+n - d+2 | Xl = 0}, (35) 

and we see that none of the columns in Jl have nonzero element in the top 
coordinate. This proves the lemma. □ 
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Proof of Proposition^ (in the real coefficient case, continued). By the as- 
sumptions, we have at least one column vector in Jr with nonzero coordinate 
on the top row. By adding such a column vector to the basis of the subspace 
(|5"5|) that are generated as in Lemma [TJ we have a basis of R/™+"- d + 2 . This 
implies rank( J q (x)) — m + n — d + 2, which proves the proposition in the real 
coefficient case. 

4-2.2. The Complex Coefficient Case 

Let X* = (f m ,l, ■ ■ ■ , /0,l)9n,l) • • • ) 00,1 j /m,2) • • • j /o,2,3n,2, ■ • ■ , #0,2, 0>n-d,li • ' • > a 0,l j 

6 m _d,i, . . . , b ,i, a n -d,2, • - • , ao,2, b m -d,2, ■ • • , &o,2) with its polynomial represen- 
tation expressed as in (|14[) (note that this assumption permits the polynomials 
-F(:e) and G(x) to be relatively prime in general). To verify our claim, we show 
that we have rank( J q (x*)) = 2(m + n — d + 1) + 1 as in ([25]) . with J q (x*) as in 
(|3"i|) . Let us express J g (x*) — (Jl | Jr), where Jl and Jr are column blocks 
expressed as 

[0 \ /2 •'-!;! 2 4 -i; 2 \ 

Jl = Mi -^2 , Jr = A^i -JV a , 
\v4 2 A x ) \ N 2 Ni J 

respectively. Then, we have the following lemma. 
Lemma 2. We have rank(JL) = 2(m + n — d + 1). 

Proof. For Ai = [C m (^4)i C n (B)i], let C m (A)i be the right m — d columns of 
C m (j4)i and C n (B)i be the right n — d columns of C n (B)±. Then, we see that 
the bottom m + n — 2d rows of the matrix C = [C m (A)i C n (B)i] is equal to 
the matrix consisting of the real part of the elements of N(A, B), the Sylvester 
matrix of A(x) and B{x). By the assumption, polynomials A(x) and B{x) are 
relatively prime, and there exist no nonzero elements in C except for the bottom 
m + n — 2d rows, thus we have rank(C) = m + n — 2d. 

By the structure of C and the lower triangular structure of C' m (A)i and 
C n (B)i, we can take the left d+ 1 columns of C m (A)i or C n (B)i satisfying 
linear independence along with C, which implies that there exist a nonsingular 
square matrix T of order m + n + 2 satisfying 

A\T = R, (36) 

where R is a lower triangular matrix, thus we have rank^ij = rank(i?) = 
m + n — d + 1. 

Furthermore, by using T and i? in (1361) . we have 

u -a 2 ] ( T pi = ( i? -ir ) , (s?) 
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followed by a suitable transformation on columns on the matrix in the right- 
hand-side of (|37[) . we can make A^T to zero matrix, which implies that 




rank(J L ) = rank R -A 2 T \ = 2 • rank(i?) = 2(m + n - d + 1). 



This proves the lemma. □ 

Proof of Proposition^ (in the complex coefficient case, continued). By the 
assumptions, we have at least one nonzero coordinate in the top row in Jr, while 
we have no nonzero coordinate in the top row in Jl, thus we have rank( J q {xj) = 
2(m + n— d + 1) + 1, which proves the proposition in the complex coefficient 
case. □ 

Remark 2. Proposition [T] says that, so long as the search direction in the mini- 
mization problem satisfies that corresponding polynomials have a GCD of degree 
not exceeding d, then J q (x) has full rank, thus we can safely calculate the next 
search direction for approximate GCD. On the other hand, it is still not clear 
when J q {x) becomes singular in our minimization problem. Although our exper- 
iments have shown that the iteration converges for any d satisfying < d < n in 
many examples, its theoretical property deserves further investigation (see also 
concluding remar ks (Section EJ). 

4-3. Setting the Initial Values 

At the beginning of iterations, we give the initial value Xq by using the 
singular value decomposition (SVD) (0), as follows. 

If.. 3.1. The Real Coefficient Case 

In the case of the real coefficients, we calculate the SVD of the (d — 1)- 
th subresultant matrix N d ^(F,G) : R™+«- 2d + 2 _> R ™+»-d+i ( see (jgj) Let 
N d -i(F, G) = U E l V be the SVD of N d -x(F, G), where 



N d ^(F, G) = UE*V, U= . . . , u m+n _ 2d+2 ), 

E = diag(<7i, . . . , C m+rl _2d+2), V = (vi, . . . ,V m+n -2d+2), 



(38) 



with Uj G R m +n - d+1 , Vj G R™+«- 2d + 2 , and E = diagfo, . . . , a m+ „_ 2d+2 ) 
denotes the diagonal matrix whose the j-th diagonal element is dj. Note that 
U and V are orthogonal matrices. Then, by a property of the SVD (0, Theo- 
rem 3.3]), the smallest singular value cr m + n —2d+2 gives the minimum distance 
of the image of the unit sphere g(" l +«- 2d + 2 )- 1 ; given as 

s (m+n-2d+2)-l = | x £ R m+n-2d+2 | = ^ 

by Nd-i, represented as 

AT ofm+n — 2d+2) — 1 r t\t | T>m+ri— 2d+2 ||„|| 1 ~l 

JVd-i ' S T ; = {V d _ia; x e R , ||x|| 2 = 1}, 
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from the origin, along with o- m + n - 2d + 2 u m+n ^ 2d+2 as its coordinates. By 
we have 

N d -i ■ V m+n -2d+2 — 0~m+n-2d+2U m +n—2d+2, 

thus Vm+n— 24+2 represents the coefficients of A(x) and B(x): let 

V m +n-2d+2 d> ■ ■ - > bo), 

A(x) = a n ^ d x n - d + ---+a Q x°, 
B{x)=b m - d x m - d + --- + boX°. 

Then, A(x) andB(x) give the least norm of AF+BG satisfying = 1 

by putting A(x) = A(x) and B{x) — B(x). 

Therefore, we admit the coefficients of F, G, A and B as the initial values 
of the iterations as 

, go, a n -d, ■ ■ ■ , ao, b m -d, ...,bo). (39) 
4-3.2. The Complex Coefficient Case 

( N -No 

In the complex case, we calculate the SVD of N = I Ar 

^JV 2 iVi 

N=UE t V, U= (u 1? . . . ,u 2(m+n _ 2d+2) ), ^ 
£ = diag(o-i,...,cr 2 ( ro+n _2d+2)) ! V — (v 1: . . . , v 2(m+n _ 2d+2 )), 

where Uj £ B, 2( - m+n - d+1 \ v . e R 2(m+n-2d+2) ) and u and y are orthogonal 
matrices. Then, as in the case of the real coefficients, the smallest singular 
value cr 2 ( m+ „_2d+2) gives the minimum distance of the image of the unit sphere 

s 2(m+„-2d+2)- 1) giyen ag 

s 2(m+n-2d+2 } -l = {j , g R 2( m +n-2d+2) | = ^ 

by N, represented as 

N • g2(m+n-M+2)-i = {Nx | x e R 2 (™+™- 2d + 2 ), || x || a = l}, 

from the origin, along with cr 2 ^ m+n _ 2d+2) u 2{m+n _ 2d+2) as its coordinates. By 
PT]|) . we have 

J* ' v 2(m+n-2d+2) = °~2(m+n-2d+2) u 2(m+n-2d+2) , 

thus «2(m+n-2d+2) represents the coefficients of A(x) and B(x): let 

t, 2(m+n-2d+2) = • • • j ^0,1) b m -d,l, ■ ■ ■ , bo,l,a n ~d,2, ■ ■ ■ , So, 2, bm-d,2, ■ ■ ■ , bo,2), 

= (a„_d,i + a n _ dt2 i)x n ~ d H h (a ,i + a , 2 i)a; , 

B(jb) = (bm-d,! + b m - d . 2 i)x m - d + ■■■ + (6 ,i + &o,2*)2; . 

Then, A(x) and B(x) give the least norm of AF+BG satisfying = 1 

by putting A(x) = A(x) and B(x) = B(x) in (Til)) , 
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Therefore, we admit the coefficients of F, G, A and B as the initial values 
of the iterations as 

x = (/m,l) ' ■ • ' fo,l,9n,l, ■ ■ ■ ,go.l:fm,2, ■ ■ ■ , /o,2, <?n,2, • • ■ , 90,2, 

0> n -d,l, ■ ■ ■ , ^0,1, b m -d,l, • • • j &0,lj On-ti,2, • • ■ , Go, 2, b m _d,2, ■ ■ ■ , &0,2j- (41) 

Regarding the Minimization Problem as the Minimum Distance (Least Squares) 
Problem 

Since we have the object function / as in (fT2")l or (f2"31 in the case of the real 
or the complex coefficients, respectively, we have V/(a;) = 2vr, where 

= t (Xl - fm, • • • ) - /o, £m+2 - 5n, • • • > X m+n +2 - .90, 0, . . . , 0), (42) 

in the case of the real coefficients, or V/(x) = 2vq, where 

V C = {x\ — fm,l, ■ ■ ■ , %m+l — /o,l) %m+2 — 9n.,lj • ■ • , ^m+n.+2 _ .90,1' 

— /m,2j • ■ • , ^2m+n+3 — J0,2) 
^2m+n+4 — .9n,2, • ■ • , ^2(m+n+2) — ff0,2, 0, . . . , 0), (43) 

in the case of the complex coefficients, respectively. However, we can regard 
our problem as finding a point x £ V q which has the minimum distance to 
the initial point Xq with respect to the (xi, . . . , x m+n+ 2)-coordinates in the 
case of the real coefficients or the {x\, . . . , X2(m+n+2) ^coordinates in the case 
of the complex coefficients, respectively, which correspond to the coefficients in 
F(x) and G(x). Therefore, in the gradient projection method at a; € V ? , the 
projection of — V/(as) in ([2"5]l should be the projection of in the case of the 
real coefficients, or vq in the case of the complex coefficients, respectively, onto 
T x , where dr and vq are as in (|4*2"j) and (|4*3"|) . respectively. These changes are 
equivalent to changing the objective function as f(x) = \f{x) then solving the 
minimization problem of f(x), subject to q(x) = 0. 

4-5. Calculating the Actual GCD and Correcting the Deformed Polynomials 

After successful end of the iterations in Algorithms [1] or [21 we obtain the 
coefficients of F(x), G(x), A(x) and B(x) satisfying (0| with A(x) and B(x) 
are relatively prime. Then, we need to compute the actual GCD H(x) of F(x) 
and G(x). Although H can be calculated as the quotient of F divided by B or 
G divided by A, naive polynomial division may cause numerical errors in the 
coefficient. Thus, we calculate the coefficients of H by the so-called least squares 
division ([3l[), followed by correcting the coefficients in F and G by using the 
calculated H, as follows. 

4- 5.1. Calculating Candidates for the GCD in the Real Coefficient Case 

For polynomials F, G, A and B represented as in ([5]) and H represented as 

H{x) = h d x d H h h x°, 
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solve the equations HB = F and HA = G with respect to H as solving the 
least squares problems of linear systems 

C d (A) t (h d ,...,h )= t (g n ,...,go), (44) 
C d (B) t (h d ,...,h ) = t (f m ,...J ), (45) 

respectively. Let H\{x),H 2 (x) £ R[x] be the candidates for the GCD whose 
coefficients are calculated as the least squares solutions of (|4"4")) and (l4"5j) , respec- 
tively. 

4-5.2. Calculating Candidates for the GCD in the Complex Coefficient Case 
For polynomials F, G, A and B represented as in (TT4|) and H represented as 

H(x) = (h d ,i + h dt2 i)x d H h (/i ,i + h 0t2 i)x°, 

solve the equations HB = F and HA = G with respect to H as solving the 
least squares problems of linear systems 

C d {AY{h d ^ +h dt 2i,---,h 0> i +/lo,2«) = *(.9n,l +ffn,2*, • • • ,ffo,i +30,2*), (46) 
C d {B) t (h d ^ 1 +h dt2 i,...,h 0il +h , 2 i) = + /™,2*, ■ • • , /o,i + /o,2«), (47) 

respectively. Then, we transfer the linear systems (|46|) and (j47|) . as follows. For 
(j47|) . let us express the matrices and vectors as the sum of the real and the 
imaginary part of which, respectively, as 

C d {B) =B X + iB 2 , 
t {h d ,i + h d ,i%, Vi + =hi + ih 2 , 

*(/m,l + /m,2», ■ • ■ , /o,l + /0,2i) = /l + ih- 

Then, (j45|) is expressed as 

(B 1 + iB 2 )(h 1 + ih 2 ) = (fi + if2)- (48) 
By equating the real and the imaginary parts in Eq. (1481) , respectively, we have 
{B l h l -B 2 h 2 ) = f x , {B x h 2 + B 2 hx) = f 2 , 



or 




Thus, we can calculate the coefficients of H(x) by solving the real least squares 
problem (|49")) . We can solve (|46| similarly. Let Hi(x), H 2 (x) e C[x] be the 
candidates for the GCD whose coefficients are calculated as the least squares 
solutions of (|4"f?|) and (|47p . respectively. 
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4-5.3. Choosing the GCD and Calculating the Deformed Polynomials 

Let Hi(x), H 2 {x) £ C[x] be the candidates for the GCD calculated as in the 
above. Then, for i — 1,2, calculate the norms of the residues as 

n = WF-HiBWl + WG-HiAWl 

respectively, and set the GCD H(x) be Hi{x) giving the minimum value of i\. 
Finally, for the chosen H(x), correct the coefficients of F(x) and G(x) as 

F(x) = H(x) ■ B{x), G{x) = H{x) ■ A{x), 

respectively. 

4-6. The Algorithm 

Summarizing the above, the algorithm for calculating approximate GCD 
becomes as follows. 

Algorithm 3 (GPGCD: Approximate GCD by the Gradient-Projection Method). 

• Inputs: 

— F(x),G(x) £ R[x] or C[x] with deg(F) > deg(G) > 0, 

— d £ N: the degree of approximate GCD with d < dcg(G), 

— e > 0: a threshold for terminating iteration in the gradient-projection 
method, 

— u £ N: an upper bound for the number of iterations permitted in 
the gradient-projection method. 

• Outputs: F(x), G(x), H{x) £ R[x] or C[x] such that F and G are deforma- 
tions of F and G, respectively, whose GCD is equal to H with deg(H) = d. 

Step 1 [Setting the initial values] As the discussions in Section 1431 set the 
initial values Xq as in (|39[) in the case of the real coefficients, or ([41]) in 
the case of the complex coefficients, respectively. 

Step 2 [Iteration] As the discussions in Section 14.41 solve the minimization 
problem of f(x) = |/(a;), subject to q(x) = 0, with f(x) and q(x) as in 
(TT2"j) and (TIB")) in the case of the real coefficients, or in (j2"3"]l and (jM)) in the 
case of the complex coefficients, respectively. Apply Algorithm Q] or [2] for 
the minimization: repeat iterations until the search direction dk (as in (|28|) 
in the gradient-projection method or in (|31j) in a modified Newton method, 
respectively) satisfies ||dfe||2 < e, or the number of iteration reaches its 
upper bound u. 

Step 3 [Construction of F, G and H] As the discussions in Section 14.51 
construct the GCD H(x) and correct the coefficients of F(x) and G{x). 
Then, return F(x), G(x) and H(x). If Step 2 did not end with the number 
of iterations less than u, report it to the user. 
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4-7. Preserving Monicity 

While Algorithm [3] permits changing the leading coefficients for calculating 
F(x) and G(x), we can also give an algorithm restricting inputs F(x) and G(x) 
and outputs F(x) and G{x) to be monic as follows. 

4-7. 1. The Real Coefficient Case 

Let F(x) and G(x) be represented as in (0 with f m = g n = 1, then, by Eq. 
(JTJ), we have b m -d = —a-n-d- Thus, we eliminate the variables f m , g n and b m -d> 
which cause the following changes. 

Changes on the Subresultant Matrix. By eliminating the variables as in the 
above, we see that Eq. ([7]) is equivalent to 



N'd-ii^i ' t ( a n-d, ■■■,ao, b m -d-i, ■ ■ ■ , bo) 
where A f ^_ 1 (F, G) is defined as 



= 0, 



K^(F,G) = 



m — 1 



fn 



1 

9n-l 



fo — 9n-r 



1 



fa 



V 



/m-l 90 



fo 



■ 1 

9n-l 

go 



J 



with (in the first column) gj = for j < 0, by subtracting the first column by 
the (n — d + l)-th column, then deleting the first row and the (n — d + l)-th 
column (corresponding to the & m _d term) in Nd-i(F,G). 

Changes on the Settings in the Minimization Problem. In solving the minimiza- 
tion problem, we substitute the variables 

(/m-l) • • - j /(J) <?n-l> • • • j£to) Q>n-di ■ ■ ■ , a 0, &m-d-l) • • • , &o) 

as x = (xi, . . . , X2(m+n~d)+i), instead of (|11J>. As a consequence, in contrast to 
(JT5J), the objective function /(as) becomes as 



/(*) - ( Xl - f m ^f +--- + (x m - /o) 2 

+ (Xm+1 ~ 9n-l) 2 H K {Xm+n ~ 9of 

Also, in contrast to © and (|TU1) . the constraints g(a?) become as 

90 = 2a 2 _ d + a\_ d _ x • • ■ + a§ + + • ■ ■ + &o - 1 = °> 

91 = (/m-l — <?n-l) a n-d + a n-d-l + &m-d-l = 0, 



(50) 



(51) 
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Changes on the Initial Values. Let N' d _ 1 = U^V be the SVD of N^^F, G), 
with 



V = («i, . . .,« m+r , 



-2d-l 



), 



"m+n-2d-l — (a>n-di • ■ • 7 a 0, &m-d-l) • • • j fy))- 



Then, in contrast to (f3T))) . the initial values become as 
x o = (/m— 1) • • • j Jo j 



,60). (52) 



TTie Algorithm. Summarizing discussions in the above, for preserving -F(x) and 
G(a;) to be monic, we modify Algorithm [3] as follows. 

Algorithm 4 (GPGCD preserving monicity, with real coefficients). Change 
Steps 1 and 2 in Algorithm [3] as follows. 

Step 1 [Setting the initial values] Set the initial values Xq as in (|52"1) . 

Step 2 [Iteration] Solve the minimization problem of f(x) = h(x), subject 
to q(x) = 0, with f(x) and q(x) defined as in ([50)1 and (l5Tj) . respectively, 
as Step 2 in Algorithm [3l 

4-7.2. The Complex Coefficient Case 

Let and G(x) be represented as in (IT41 with / m ,i = = 1 and 

/m,2 = 9n,2 = 0, then, by Eq. 0, we have & m _dj = -a n -d,j for J € {1,2}. 
Thus, for j € {1,2}, we eliminate the variables f m ,j, 9n,j and b m -dj, which 
cause the following changes. 

Changes on the Subresultant Matrix. By eliminating the variables as in the 
above, we see that Eq. (fT6|) is equivalent to 



/ (/m-1,1 — 571-1,1) + (/m-1,2 — ffn-1,2)* 1 

: /m-1,1 + /m-1,2* 

(/o,l — 5n-m,l) + (/o,2 _ <7n-m,2)* '■ 

fo,l + /o,2* 



1 

/m-1,1 + /m-1,2* 



071-1,1 + 071-1,2' 



,1 + .90,2* 



1 

ffn-1,1 + Sn-1,2* 



/o,l + /o,2* 



oo,i + O0,2* 
&m-rf-l,l + &m-d-l,2* 



.90,1+30,2* / 



\ &0,1 + &0,2* / 

-0, (53) 
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with (in the first column of the matrix in the left-hand-side) g^j — for i < and 
j € {1,2}, in which the matrix in the left-hand-side is obtained by subtracting 
the first column by the (n — d + l)-th column, then deleting the first row and 
the (n — d + l)-th column (corresponding to the &m-d,i + frm-d,2i term) in the 
corresponding matrix in (|16p . Then, Eq. (I17p becomes as 



with 



(N[+N^i)(v' 1 +v , 2 i) = 0, 



f /m-1,1 — ffn-1,1 1 

; /m-1,1 
/o,l — <?n-m,l '■ 



1 : •• 1 

/m-1,1 30,1 ffn-1,1 



V 



/«». 



^2 



/ /m-1,2 — ffn-1,2 1 

: /m-1,2 
/o,2 — <?n-m,2 '■ 



1 

?n-l,2 



V 

t>i = '(a„_d,i, 

«2 = *(an-rf,2, 



/o,2 



, £&0,lj &m-d-l,lj 
, «0,2, frm-d-1,2, 



/o,2 
■ - ,60,1), 
•,&0,2)- 



3o,i / 
\ 



/m-1,2 30,2 ffn-1,2 



.90,2 



(54) 



Changes on the Settings in the Minimization Problem. In solving the minimiza- 
tion problem, we substitute the variables 



(/m-1,1, • • • , /o,l, <?n-l,l, • • • ) #0,1) /m-1,2, • • • , /o,2, ffn-1,2, • • • , 50,2, 
O-n-d,!, ■ • ■ , ^0,1, Om-d-1,1) ■ • • , &0,1, 0-n~d,2, • ■ • , Go, 2, frm-ci-1,2, 



■ ,&0,2) 



as a; = (a?i, . . . , X4(m+n-d+i)), instead of (|22|) . As a consequence, in contrast to 
|), the objective function /(a;) becomes as 



f(x) ={ Xl -/m-l,l) 2 + '-- + (a 



/o,i) 2 



+ (z m +l - ffn-l.l) 2 H h (^m+n - .90, l) 2 

+ (^m+n+l — /m-1,2) 2 + ' ' ' + (£2m+n _ /o,2) 2 
+ (l2m+n+l - ffn-1,2) 2 H + (^2(m+n) ~ 50,2) 2 - 

The constraints becomes as follows. Now, Eq. (fl"9]l becomes as 



(55) 



(56) 
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with N[, N%, v[ and v' 2 are defined as in Furthermore, the constraint for 

the coefficients in A(x) and B(x) as in ([20)) now becomes as 

\\A(x)\\l + = (2^_ d>1 + • • • + a^) + (b 2 m _ d _ hl + ■■■+ b\ x ) 

+ (K-d,2 + ■■■ + «o, 2 ) + (bl- d -i,2 + ■■■ + b 2 , 2 ) - 1 - 0. (57) 

Then, by the same way we have constructed (f2"Tj) , we put ([55)) and ([57)) together 
as 
























and we obtain the constraint q(a;) = as 

g(x) = t {qi(x),... : q 2{m+ 
where qj (x) corresponds to the j-th row of matrix- vector product in ([58 







A?)-* 


(58) 






-d)+iN) = o, 


(59) 



/AT' —AT' 

Changes on the Initial Values. Let AT' = U S * F be the SVD of N' = I ^ 
with 

V = (t?l, . . . , t>2(m+n-2<i+l)) I 



v 2(m+n-2d+l) — (fln-rf.lj • • • j a 0,l) &m-d-l,li • • • j 



0>n-d,2i • ■ ■ , Q.0,2, &ro-d-l,2j ■ ■ • , &0,2j- 

Then, in contrast to (14D). the initial value becomes as 



x — (/m-l,lj ■ ■ • i /o,lj <7n-l,lj ■ ■ • i ff0,lj /m-1,2, • • • j /o,2, <7ra-l,2, • ■ ■ , 30,2, 
a n -d,l! ■ ■ • j ao : lj bm-d-1,1, ■ ■ • , &0,lj »n-<12, ■ ■ • , 00,2, &m— d-1,2; ■ ■ • , &0,2)- (60) 

XTie Algorithm. Summarizing discussions in the above, for preserving F(x) and 
G(x) to be monic, we modify Algorithm [3] as follows. 

Algorithm 5 (GPGCD preserving monicity, with complex coefficients). Change 
Steps 1 and 2 in Algorithm [3] as follows. 



Step 1 [Setting the initial values] Set the initial values Xq as in ([60)) . 

Step 2 [Iteration] Solve the minimization problem of f(x) = subject 
to q(x) = 0, with f(x) and q(x) defined as in ([55)) and (|59)) . respectively, 
as Step 2 in Algorithm [3) 

4-7.3. Running Time Analysis 

We give an analysis for running time of Algorithm [3] with employing the 
modified Newton method. 

In Step 1, we set the initial values by the SVD. Since the dimension of 
subresultant matrix is 0(m + n — d), running time in this step becomes 0((m + 
n-d) 3 ). 
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In Step 2, we estimate running time just for one iteration, since the number of 
iterations for convergence of solution may vary depending on the given problem. 
This step essentially depends on solving the linear system (f5U|) with the Jacobian 
matrix J q (xk) defined as in (1521 for the real coefficient case or in for the 
complex coefficient case. In both cases, dimension of J q {xk) is 0(m + n — d), 
thus we can estimate running time for solving the linear system as 0((m + 
n-d) 3 ). 

Step 3 depends on calculating the least square solution of the linear system 
as in Section l4~5l whose running time becomes as 0{{m + n — d) 3 ). 

As a consequence, we can estimate running time of Algorithm [3] as 0((m + 
n — d) 3 ) times the number of iterations for finding a GCD. 

4-8. Examples 

Now we show examples of Algorithm [3] in the case of the real coefficients 
(more comprehensive experiments are presented in the next section). 

Note that, for the minimization method, we have employed a modified New- 
ton method (Algorithm [2]) . Computations in Example [T] have been executed on 
a computer algebra system Mathematica 6 with hardware floating-point arith- 
metic, while those in Examples [5] and [3] have been executed on another computer 
algebra system Maple 15 with Digits=10. 



Example 1. This example is given by Karmarkar and Lakshman [15j, followed 
by Kaltofen et al. [14]. Let F (x),G(x) <E R[x] be 

F(x) = x 2 - 6x + 5 = (x - l)(x - 5), 

G(x) = x 2 - 6.3a; + 5.72 = (a; - 1.1) (a; - 5.2), 

and find F(x), G(x) £ R[x] which have the GCD of degree 1, namely F(x) and 
G(x) have one common zero. 

Case 1: The leading coefficient can be perturbed. Applying Algorithm [3J 
to F and G, with d — 1 and e = 1.0 x 10~ 8 , after 7 iterations, we obtain the 
polynomials F and G as 

F(x) = 0.985006cc 2 - 6.00294x + 4.99942, 
G(x) = 1.01495x 2 - 6.29707x + 5.72058, 



with perturbations as y\\F — F\\ 2 + \\G — G\\ 2 = 0.0215941 and the common 

zero of F(x) and G(x) as i = 5.09 890419203. In Kaltofen et al. [3, the 
calculated perturbations obtained is a/0. 0004663 = 0.021594 with the common 
zero as x = 5.09890429. Karmarkar and Lakshman }15j only give an example 
without perturbations on the leading coefficients. 

Case 2: The leading coefficient cannot be perturbed. Applying Algorithm[3] 
(preserving monicity) with the same arguments as in Case 1, after 7 iterations, 
we obtain the polynomials F and G as 

F(x) = x 2 - 6.07504x + 4.98528, 
G(x) = x 2 - 6.22218a; + 5.73527, 
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with perturbations as y\\F — F\\?,+ ||G — G||| = 0.110164 And the common 

zero of F(x) and G(x) as x = 5.0 969464650. In Kaltofen et al. 0, the calcu- 
lated perturbations obtained is yO- 01213604583 = 0.110164 with the common 
zero as x = 5.0969478. In Karmarkar and Lakshman [IB] , the calculated per- 
turbations obtained is \/0-01213605293 = 0.110164 with the common zero as 
x = 5.096939087. 

The next examples, originally by Sanuki and Sasaki (2lj . are ill-conditioned 
ones with the small or large leading coefficient GCD. 



Example 2 (A small leading coefficient problem 21, Example 4]). Let F{x) 
and G{x) be 

F(x) = (x i + x 2 +x + l)(0.001a; 2 + x + 1), 
G(.x) = (a; 3 + x 2 + x + l)(0.001a; 2 +X + 1). 

Applying Algorithm [3] to F and G, with d = 2 and e = 1.0 x 10~ 8 , after 1 
iteration, we obtain the polynomials F, G and H as 

F(x) ~ G(a:) ~ G(^), 

= O.OOlx 2 + 0.99999999362; + 0.9999999936, 



with \ \\F - + ||G - G||| = 8.485281374 x 10" 



Example 3 (A large leading coefficient problem [211 . Example 5]). Let F(x) 
and G(x) be 

F(x) = (a; 6 - 0.00001(0.8x 5 + 3x 4 - 4a; 3 - 4a; 2 - 5a; + 1)) • C[x), 
G(x) = (x 5 +x 4 + x 3 - 0.1a; 2 + 1) • C(x), 

with C(x) = x 2 + 0.001. Applying Algorithm [3] to F and G, with d = 2 and 
£ = 1.0 x 10~ 8 , after 1 iteration, we obtain the polynomials F, G and H as 

F(x) ~ G(a?) ~ G(»), 

= a; 2 + 1.548794164 x 10 _16 a; + 0.001, 



with A/llf-Fll^ + IIG-GII 2 = 1.735004369 x 10" 14 . 



5. Experiments 

We have implemented our GPGCD method (Algorithm |3[ on a computer 
algebra system Mapl^l and carried out the following tests: 



The implementation is available at Project Hosting on Google Code [27 



25 



1. fSection l5.1[) Comparison of performance of the gradient-projection method 
(Algorithm [TJ and the modified Newton method (Algorithm [2]) on ran- 
domly generated polynomials with approximate GCD, 

2. (Section I5.2j) Comparison of performance of the GPGCD method with 
a method based on the structured total least norm (STLN) method by 
Kaltofen et al. [l3[ and the UVGCD method by Zeng [3l| on large sets of 
randomly-generated polynomials with approximate GCD, 

3. (Section |5 . 3[) Comparison of performance of the GPGCD method with the 
STLN-based method and the UVGCD method on ill-conditioned polyno- 
mials and other test cases by Zeng [3l[ and Bini and Boito 




Note that, in Test[U we have tested both the cases of the real and the complex 
coefficients, while, in the other tests, we have tested only the case of the real 
coefficients. 

In Tests Q] and [H we have generated random polynomials with GCD then 
added noise, as follows. First, we have generated a pair of monic polynomials 
Fo(x) and Gq(x) of degrees m and n, respectively, with the GCD of degree d. 
The GCD and the prime parts of degrees m — d and n — d are generated as 
monic polynomials and with random coefficients c e [—10, 10] of floating-point 
numbers. For noise, we have generated a pair of polynomials F^(x) and Gn(i) 
of degrees m — 1 and n — 1, respectively, with random coefficients as the same 
as for Fo(x) and G$(x). Then, we have defined a pair of test polynomials F(x) 
and G{x) as 



respectively, scaling the noise such that the 2-norm of the noise for F and G 
is equal to eF and ec, respectively. In the present test, we set eF = &g = 0.1. 
(See also the notes in Section 1531 ) 

The tests have been carried out on Intel Core2 Duo Mobile Processor T7400 
(in Apple MacBook "Mid-2007" model) at 2.16 GHz with RAM 2GB, under Mac 
OS X 10.6. All the tests have been carried out on Maple 15 with Digits=15 
executing hardware floating-point arithmetic. 

5.1. Test\j\' Comparison of the Gradient- Projection Method and the Modified 
Newton Method 

In this test, we have compared performance of the gradient-projection method 
(Algorithm [T|) and a modified Newton method (Algorithm [2]) , only in the case 
of the real coefficients. For every example, we have generated one random 
test polynomial as in the above, and we have applied Algorithm [3] (preserving 
monicity) with u = 100 and e = 1.0 x 10~ 8 . 

Table Q] shows the result of the test: m and n denotes the degree of a tested 
pair F and G, respectively, and d denotes the degree of approximate GCD; 
"Perturbation" is the perturbation of the perturbed polynomials from the initial 
inputs, calculated as 




F(x) =F {x) + 



\\Mx)h 



F N (x), G(x)=G (x) + 



\\Gn(x)\\ 2 



Gn(x), 




(61) 
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Ex. 


m, n 


d 


Perturbation 
(EJ) 


^Iterations 


Time (sec.) 


Alg.lU 


Alg.|2j 


Alg.lU 


Alg.[2| 


1 


10,10 


5 


4.25e-2 


3 


4 


0.10 


0.04 


2 


20,20 


10 


6.86e-2 


3 


4 


0.17 


0.11 


3 


40,40 


20 


6.80e-2 


4 


5 


0.61 


0.16 


4 


60,60 


30 


7.24e-2 


3 


4 


0.69 


0.23 


5 


80,80 


40 


5.06e-2 


3 


4 


1.41 


0.41 


6 


100, 100 


50 


7.26e-2 


3 


4 


2.21 


0.76 



Table 1: Test results comparing the gradient-projection method and the modified Newton 
method; see Section [571] for details. 



where "aeb" with a and b as numbers denotes a x 10 b ; "^Iterations" is the 
number of iterations; "Time" is computing time in seconds. The columns with 
"Alg. U' and "Alg. [2]' are the data for Algorithm Q] (the gradient-projection 
method) and Algorithm [2] (the modified Newton method), respectively. Note 
that, the "Perturbation" is a single column since both algorithms give almost 
the same values in each examples. 

We see that, in all the test cases, the number of iterations of the gradient- 
projection method (Algorithm Q} is equal to 3 or 4, which is smaller than that of 
the modified Newton method (Algorithm [5]) which is equal to 4 or 5. However, 
an iteration in Algorithm [1] includes solving a linear system at least twice: once 
in the projection step (Step 2) and at least once in the restoration step (Step 
3); whereas an iteration in Algorithm [5] includes that only once. Thus, total 
number of solving a linear system in Algorithm [2] is about a half of that in 
Algorithm Q] Furthermore, computing time shows that the modified Newton 
method runs approximately twice as fast as the gradient projection method. 
Therefore, we adopt Algorithm [5] as the method of minimization in the GPGCD 
method (Algorithm^. 

5.2. Test\^ Tests on Large Sets of Randomly-generated Polynomials 

In this test, we have compared Algorithm [3]with a method based on the 
structured total least norm (STLN) method ([13]) and the UVGCD method 
([13]): using their implementation for the Maple, in the both cases of the 
real and the complex coefficients. In our implementation of Algorithm [3J we 
have chosen the modified Newton method (Algorithm [5]) for minimization. In 
the STLN-based method, we have used their procedure R_con_mulpoly and 
C_con_mulpoly, which calculates the approximate GCD of several polynomials 
in R[x] and C[i], respectively. In the UVGCD method, we have used their 
procedure uvgcd for calculating approximate GCD of polynomials in R[x] and 
C[x}. 

Note that, in this test, we have defined test polynomials satisfying another 
requirement: to make sure that the input polynomials F(x) and G(x) do not 
have a GCD of degree exceeding d, we have adopted only those satisfying that 
the smallest singular value of the d-th subresultant matrix N d (F, G) (see (0)) is 
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Ex. 


Tfl 71 


d 


Perturbation 




Time (sec.) 


^Iterations 


STLN 


UVGCD 


GPGCD 


STLN 


UVGCD 


GPGCD 


STLN 


GPGCD 


1 


10,10 


5 


5.64e-2 


1.79e-l 


5.64e-2 


0.38 


0.64 


0.04 


4.46 


4.50 


2 


20,20 


10 


6.22e-2 


1.85e-l 


6.22e-2 


1.16 


0.86 


0.06 


4.40 


4.40 


3 


30,30 


15 


6.65e-2 


1.87e-l 


6.65e-2 


2.43 


1.34 


0.10 


4.37 


4.46 


4 


40,40 


20 


6.48e-2 


1.96e-l 


6.48e-2 


4.05 


2.09 


0.13 


4.11 


4.15 


5 


50,50 


25 


6.91e-2 


1.91e-l 


6.91e-2 


6.30 


3.34 


0.19 


4.03 


4.16 


6 


60,60 


30 


6.75e-2 


1.94e-l 


6.75e-2 


9.09 


4.37 


0.26 


4.00 


4.18 


7 


70,70 


35 


6.89e-2 


2.08e-l 


6.89e-2 


12.47 


5.71 


0.35 


3.96 


4.13 


8 


80,80 


40 


6.78e-2 


1.91e-l 


6.78e-2 


16.95 


7.95 


0.44 


3.16 


4.11 


9 


90,90 


45 


6.92e-2 


1.95e-l 


6.92e-2 


22.09 


10.20 


0.57 


3.96 


4.10 


10 


100, 100 


50 


6.98e-2 


1.95e-l 


6.98e-2 


27.48 


13.02 


0.69 


3.88 


4.09 



Table 2: Test results for large sets of polynomials with approximate GCD, in the case of the 
real coefficients; see Section 15.21 for details. 



larger than or equal to lH 

For every example, we have generated 100 random test polynomials as in the 
above. In executing Algorithm [31 we have set u — 200 and e = 1.0 x 10~ 8 ; in 
R_con_mulpoly and C_con_mulpoly, we have set the tolerance e = 1.0 x 10~ 8 ; 
in uvged, we have set the initial tolerance S = 1.0 X 10~ 2 and have changed it 
until we have obtained an approximate GCD of desired degree. 

Tables [2] and [3J show the results of the test in the case of the real and the 
complex coefficients, respectively: m and n denotes the degree of a pair F and 
G, respectively, and d denotes the degree of approximate GCD. The columns 
with "STLN" are the data for the STLN-based method; "UVGCD" are the 
data for the UVGCD method; "GPGCD" are the data for the GPGCD method 
(Algorithm [3]) . "Perturbation" , "^Iterations" and "Time" are the same as 
those in Table [TJ respectively. (Note that computing time for the UVGCD 
method does not include the time for "try and error" calculations by changing 
the tolerance 5: it is just for successful calculations.) 

We see that the average of magnitude of perturbations by the GPGCD 
method is as small as that by the STLN-based method, which is approximately 
one-tenth as large as that by the UVGCD method. For computing time, the 
GPGCD method calculates approximate GCD very efficiently, faster than the 
STLN-based method by approximately from 10 to 30 times and the UVGCD 
method by approximately from 6 to 10 times. 

Remark 3. In this experiment, we have compared our implementation designed 



3 Our previous test results ( |26| . Section 5.2]) have shown that there were test cases (input 
polynomials with the real coefficients) in which the GPGCD method was not able to calcu- 
late an approximate GCD with sufficiently small magnitude of perturbations. After thorough 
investigation, we have found that such input polynomials accidentally have an approximate 
GCD of degree exceeding d. Thus, in the present test, we have denned totally new test poly- 
nomials satisfying the above requirement, then none of such phenomena have been observed 
with the test. 
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Ex. 


Tfl 71 


d 


Perturbation 




Time (sec.) 


^Iterations 


STLN 


UVGCD 


GPGCD 


STLN 


UVGCD 


GPGCD 


STLN 


GPGCD 


1 


10,10 


5 


5.92e-2 


1.54e-l 


5.92e-2 


1.58 


0.64 


0.11 


4.50 


4.46 


2 


20,20 


10 


6.40e-2 


1.41e-l 


6.40e-2 


5.34 


1.31 


0.20 


4.30 


4.30 


3 


30,30 


15 


6.63e-2 


1.40e-l 


6.63e-2 


11.63 


2.12 


0.35 


4.21 


4.24 


4 


40,40 


20 


6.61e-2 


1.34e-l 


6.61e-2 


21.57 


3.51 


0.55 


4.15 


4.13 


5 


50,50 


25 


6.86e-2 


1.48e-l 


6.86e-2 


34.23 


5.03 


0.83 


4.06 


4.10 


6 


60,60 


30 


6.86e-2 


1.51e-l 


6.86e-2 


50.40 


7.39 


1.16 


4.02 


4.05 


7 


70,70 


35 


6.94e-2 


1.41e-l 


6.94e-2 


69.54 


10.31 


1.56 


3.93 


4.05 


8 


80,80 


40 


6.85e-2 


1.44e-l 


6.85e-2 


93.77 


14.01 


2.07 


3.91 


4.07 


9 


90,90 


45 


6.84e-2 


1.52e-l 


6.84e-2 


122.97 


18.30 


2.65 


3.90 


4.04 


10 


100, 100 


50 


6.94e-2 


1.65e-l 


6.94e-2 


157.02 


23.72 


3.37 


3.86 


4.04 



Table 3: Test results for large sets of polynomials with approximate GCD. in the case of the 
complex coefficients; see Section 15.21 for details. 



for problems of two univariate polynomials against the implementation of the 
STLN-based method designed for multivariate multi-polynomial problems with 
additional linear coefficient constraints. Kaltofen 12j has reported that they 
have tested their implementation for just two univariate polynomials with real 
coefficients ([3]) on an example similar to ours with degree 100 and GCD degree 
50, and it took (on a ThinkPad of 1.8 GHz with RAM 1GB) 2 iterations and 9 
seconds. This result will give the reader some idea on efficiency of our method. 



5.3. Test\^ Tests for Ill-conditioned Polynomials and Other Cases 

In this test, we have compared Algorithm [3] with the STLN-based method 
([HI) and the UVGCD method ([3(|) on some ill-conditioned polynomials and 
other test cases by Zeng |3l| and Bini and Boito Q, as follows. 

Note that we give the degree of approximate GCD in the STLN-based 
method and the GPGCD method, while we give the tolerance S then the al- 
gorithm estimates the degree of approximate GCD in the UVGCD method. 
Also note that, in some tests in this section, we have measured the relative er- 
ror of approximate GCD from the given GCD ([63]) instead of the magnitude of 
perturbation (|61[) because, in such cases, we have given test polynomials with 
predefined (approximate) GCD and have intended to observe "nearness" of the 
calculated approximate GCD from the predefined one. 

Throughout the tables in this section, the columns with "STLN" , "UVGCD" , 
"GPGCD" , "Perturbation" , and "Time" the same as those in the above, respec- 
tively. 

Example 4. An example of ill-conditioned polynomial by Zeng [3l], Test 1]. 
Let n be an even positive number and k = ro/2, and define p n — u n v n and 
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n 


Relative error of GCD (1531 


STLN 


UVGCD 


GPGCD 


6 


1.04e-14 


4.60e-15 


3.68e-15 


8 


3.98e-13 


7.90e-13 


4.30e-13 


10 


1.08e-10 


7.89e-12 


1.08e-10 


12 


2.87e-10 


2.95e-ll 


2.94e-10 


14 


3.10e-9 


3.65e-10 


3.14e-9 


16 


6.22e-9 (*1) 


3.83e-10 


8.00e-9 


18 


1.38e-6 (*1) 


9.68e-9 


1.36e-6 


20 


6.95e-6 (*1) 


1.21e-8 


7.11e-6 (*2) 



Table 4: Test results for test polynomials H62H . See Example [4] for details. 



u n w n , where 

k k 
3=1 3=1 



(62) 

[(x — T\ctj) 2 + rf/3j], ctj — cos — , (3j = sin — 

3=k+l U U 



for n = 0.5 and T2 = 1.5. The zeros of p n and q n lie on the circles of radius r% 
and T2- We had the test for n = 6, ... ,20 increased by 2. 

Table [4] shows the result of the test. "Relative error of GCD" is calculated 

by 

\\u n {x) - U n {x)\\ 2 . . 

ii 7~Tii ' ^ ' 

where w n is predefined GCD as shown in (|62p and u„ is approximate GCD. In the 
table, (*1) indicates that the STLN-based method did not converge within 50 
times of iterations which is a built-in threshold; (*2) indicates that the GPGCD 
method did not converge within 100 times of iterations. We see that, in the 
GPGCD method as well as in the STLN-based method, the number of itera- 
tions increases and the accuracy of calculated approximate GCD decreases as 
n increases. On the other hand, the UVGCD method has better accuracy of 
approximate GCD for large n. 

Example 5. Another example of ill-conditioned polynomial by Zeng |3lL Test 
2]. Let 

10 10 

p(x) = l[(x-x j ), q(x) = H(x-x j + 10-i), Xj = (-iy(j/2), (64) 
l l 

The zeros of q have decreasing distances as 0.1, 0.01, . . . , from those of p. We 
have tried to calculate an approximate GCD of degree d from 1 to 10 increased 
by 1. 
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d 


Perturbation (1611) 


STLN 


GPGCD 


1 


5.17e-l (*1) 


3.21e3 


2 


6.95e-4 (*1) 


3.06e0 


3 


1.97e-5 


1.26e0 


4 


2.89e-6 


2.25e-l 


5 


5.28e-5 


4.75e-l 


6 


2.15e-3 


2.16e-3 


7 


8.34e-2 


8.34e-2 


8 


2.04e0 


2.04e0 


9 


4.70el 


4.70el 


10 


7.73e2 


7.73e2 



Table 5: Test results for test polynomials | [64l) with the STLN-based method and the GPGCD 
method. See Example \E\ for details. 



Tables [5] and [5] show the result of the test. In this test, we have measured per- 
turbation (|61l) since p and q are pairwisely relatively prime in a rigorous sense. 
Note that we have put the results for the UVGCD method in Table El separated 
from those for the GPGCD and the STLN-based methods in Tabled because we 
have given the tolerance S to obtain approximate GCD in the UVGCD method, 
while we have given the degree d in the GPGCD and the STLN-based methods. 
In Tabled (*1) indicates that the STLN-based method did not converge within 
50 times of iterations which is a built-in threshold. 

We see that, for d > 6, all the methods find approximate GCD with simi- 
lar magnitude of perturbations. However, for smaller value of d, the UVGCD 
method finds approximate GCD with considerably smaller magnitude of pertur- 
bations than those in the other methods, followed by the STLN-based method. 

Example 6. An example with GCDs of large degree by Zeng [3l|, Test 3]. Let 

p n =u n v, q„ =u n w, 

3 3 (65) 

v(x)=^2x J , w(x) =^2(-xy , 

3=0 3=0 

where u n {x) is a GCD defined as a polynomial of degree n whose coefficients 
are random integers in the range [—5, 5] and v{x) and w(x) are fixed cofactors. 

Table [7] shows the result of the test by measuring relative error of approxi- 
mate GCD ([63| . In this test, we have also measured computing time because 
the difference of it became large among the methods for large degree of approx- 
imate GCD. We see that the UVGCD method calculates approximate GCD 
with the best accuracy, followed by the STLN-based method and the GPGCD 
method. On the other hand, the GPGCD method is more efficient than the 
other methods. 
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UVGCD 


X 





|J /"it f iivnof i /"i n t\ 


l.Oe-11 


1 


8 02e-10 

V.J . V / \ IV/ 


l.Oe-10 


2 


3.27e-8 


1.0e-9 


3 


6.03e-7 


1.0e-8 


4 


1.99e-5 


1.0e-7 


5 


3.45e-4 


1.0e-6 


5 


3.45e-4 


1.0e-5 


6 


9.61e-3 


1.0e-4 


7 


1.79e-l 


1.0e-3 


8 


3.18e0 


1.0e-2 


8 


3.18e0 


l.Oe-1 


9 


5.00el 


l.OeO 


10 


8.40e2 



Table 6: Test results for test polynomials (I64D with the UVGCD method. Sec Example [5] for 
details. 



n 


Relative error of GCD (j63l) 


Time (sec.) 


STLN 


UVGCD 


GPGCD 


STLN 


UVGCD 


GPGCD 


50 


1.60e-15 


1.04e-16 


2.63e-15 


1.77 


0.22 


0.04 


100 


1.16e-15 


1.59e-16 


4.41e-15 


8.17 


0.31 


0.06 


200 


1.14e-15 


1.06e-16 


1.23e-14 


45.09 


0.83 


0.12 


500 


1.35e-15 


1.37e-16 


1.84e-14 


552.09 


3.39 


0.64 


1000 


1.42e-15 


1.69e-16 


5.30e-14 


4318.38 


18.66 


3.27 



Table 7: Test results for test polynomials J65J. See Example [6] for details. 



Example 7. An example with multiple zeros of high multiplicities by Bini and 
Boito [2j, Example 4.5]. Let 

u k (x) = (x 3 +3x-l)(x-l) k , v k {x)=u'(x), (66) 

for positive integer k. Note that the GCD of Uk(x) and Vk(x) is Wk(x) — 

(JB-I)*- 1 . 

Table [8] shows the result of the test. In the table, as in Example HI (*1) 
indicates that the STLN-based method did not converge within 50 times of 
iterations which is a built-in threshold; (*2) indicates that the GPGCD method 
did not converge within 100 times of iterations. 

We see that, in the GPGCD method as well as in the STLN-based method, 
the number of iterations increases and the accuracy of calculated approximate 
GCD decreases for k = 35 and 45. On the other hand, the UVGCD method 
calculates approximate GCD accurately for large k. 

Example 8. Another example with multiple zeros of high multiplicities by Zeng 
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k 


Relative 


error of GCD (|55]l 




STLN 


UVGCD 


GPGCD 


15 


2.35e-13 


3.08e-15 


1.86e-12 


25 


1.64e-ll 


1.13e-14 


6.67e-ll 


35 


3.79e-10 (*1) 


8.02e-15 


3.58e-9 (*2) 


45 


4.23e-8 (*1) 


1.13e-14 


1.78e-7 (*2) 



Table 8: Test results for test polynomials (166 D - See Example [7] for details. 



[m 1 ,m2,m 3 ,m 4 ] 


Relative error of GCD ([63]) 




STLN 


UVGCD 


GPGCD 


[2,1,1,0] 


l.lle-13 


9.42e-16 


2.83e-13 


[3,2,1,0] 


7.33e-13 


3.31e-15 


8.23e-12 


[4,3,2,1] 


2.35e-9 


2.95e-13 


2.68e-9 


[5,3,2,1] 


1.89e-8 


3.38e-12 


5.56e-9 


[9,6,4,2] 


4.72e-8 (*1) 


5.31e-ll 


6.05e-8 (*2) 


[20,14,10,5] 


5.06e-l (*1) 


3.13e-10 


9.98e-l \*2) 


[80,60,40,20] 


l.OeO (*1) 


1.08e-3 


l.OeO (*2) 


[100,60,40,20] 


l.OeO (*1) 


2.16e-4 


N/A (*3) 



Table 9: Test results for test polynomials l|67ll . See Example [8] for details. 



31, Test 61. Let 



M-± M (67) 

9[m 1 ,m2,m 3) m4]i^j dx m2,m3 ' m4 ^ 

for nonnegative integers mi, . . . , m^. Note that the GCD of P[ mi ,m 2 .m 3 ,m 4 ](x) 
and g [mi , m2 , m3 , m4 ](z) is (x - l) m 'i(x - 2) m '*(x - 3) m Hx ~ A)< with = 
max{mj — 1, 0} for j = 1, . . . , 4. 

Table [9] shows the result of the test. In the table, as in Examples |4] and [3 
(*1) indicates that the STLN-based method did not converge within 50 times of 
iterations which is a built-in threshold; (*2) indicates that the GPGCD method 
did not converge within 100 times of iterations. Furthermore, (*3) indicates 
that the GPGCD method stopped abnormally because the solution of a lin- 
ear system with the coefficient matrix (the Jacobian matrix) as shown in (|32[) 
became unexpectedly large. 

We see that, in the GPGCD method as well as in the STLN-based method, 
the number of iterations increases and the accuracy of calculated approximate 
GCD becomes almost meaningless for inputs of large degree. On the other hand, 
the UVGCD method is quite stable (in the sense of convergence of the algorithm) 
and more accurate for calculating approximate GCD for those inputs. 
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6. Concluding Remarks 

We have proposed an iterative method, based on the modified Newton 
method which is a generalization of the gradient-projection method, for calcu- 
lating approximate GCD of univariate polynomials with the real or the complex 
coefficients. 

Our experiments comparing the GPGCD method with the STLN-based 
method and the UVGCD method have discovered advantages and disadvan- 
tages of these methods, as follows. In the case that input polynomials already 
have exact or approximate GCD, then the UVGCD method calculates the ap- 
proximate GCD with the best accuracy and relatively fast convergence among 
them. On the other hand, in the case that the magnitude of "noise" is larger, 
then the magnitude of perturbations calculated by the GPGCD method or the 
STLN-based method is smaller than that calculated by the UVGCD method. 
Furthermore, in such cases, the GPGCD method has shown significantly better 
performance over the other methods in its speed, by approximately up to 30 
times for the STLN-based method and 10 times for the UVGCD method, which 
seems to be sufficiently practical. Other examples have shown that the GPGCD 
method properly calculates approximate GCD with small or large leading coef- 
ficient. 

Our result have shown that, in contrast to the STLN-based methods which 
uses structure preserving feature for matrix computations, our simple method 
can achieve accurate and efficient computation as or more than theirs in cal- 
culating approximate GCDs in many examples. On the other hand, our result 
have also shown that our method is less accurate than the UVGCD method 
especially in the case the given polynomials lie sufficiently close to polynomials 
that have a GCD in a rigorous sense. These results suggest that there are some 
opportunities for improvements of accuracy and/or efficiency in calculating ap- 
proximate GCDs with optimization strategies. 

For the future research, the followings are of interest. 

• Convergence analysis of the minimizations: showing global convergence of 
local method is difficult in general (see e.g. Blum et al. [3j), as the original 
paper on the modified Newton method (;24]) only shows its stability by 
observing whether the Jacobian matrix of the constraint at a local minimal 
point has full-rank or not. However, it may be possible to analyze local 
convergence property depending on condition on the initial point and/or 
local minimal point. (See also Remarks [T] and [2J . 

• Improvements on the efficiency: time complexity of our method depends 
on the minimization, or solving a system of linear equations in each it- 
eration. Thus, analyzing the structure of matrices might improve the 
efficiency in solving a linear system. 

• Comparison with other methods (approaches) for approximate GCD: from 
various points of view such as accuracy, stability, efficiency, and so on, 
comparison of our methods with other methods will reveal advantages 
and drawbacks of our method in more detail. 
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Other topics, such as generalization of our method to several input polynomi- 
als, are also among our next problems, some of which are currently under our 
investigation (|28|). 
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